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DERIVATION OF NINTH STAGE RUNGE-KUTTA METHOD FOR THE SOLUTION OF 
FIRST ORDER DIFFERENTIAL EQUATIONS 


MSHELIA DW, BADMUS AM AND YAKUBU DG 


Abstract: We present hybrid block of eight integrators which are of uniform order nine 
through interpolation and collocation procedures. The properties of the hybrid block 
integrators are fully investigated and confirmed to be computationally stable also the block 
method derived are reconstructed to ninth stage Runge-Kutta method which implemented 
on stiff, physical and life problems. The results obtained compare favourably with the 
existing methods when implemented in Runge-Kutta mode. 


1.0 Introduction 

Among the most important mathematical tools used in producing models in the physical 
sciences, Biological sciences and Engineering are differential equations. But most of these 
differential equations do not posses closed form or finite solutions. Even if they posses 
closed form solutions we do not know the method of getting them. 

In many real-life situations, the differential equation that models the problem is too 
complicated to solve exactly. Hence there is need to develop an accurate algorithm for 
obtaining an equivalent approximating solution to the original problems. Most recent 
researchers have developed some block methods to cater for this class of problems 

у= (х,у),  ya=p (1) 
Among such researchers аге [1], [2],[3], [4] , [51461 and [7] to mention a few. 

In our method, a block of eight integrators were proposed at step length of four which are 
of uniform order nine and also all the discrete schemes in our block came from a single 
continuous formula 

Definition 1.0 Zero stable 

A linear multi-step method is said to be Zero-stable if the roots Rj, j = 1(1)k of the first 
characteristics polynomials 


k 
» АЕК"! 


ї=0 
If one of the roots is +1, we call this the principal root of p(R). [7] 


p(R) = det = 0, Ao = —1, satisfies [R;| < 1 
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2.0 Development of the method 
Our objective is to derive a block of eight integrators of the form 


Y(Xn+v) =%0 Yn + h (вол, + Baf m + ааа + Palaa + b2fn+2 + Bsf 5 + B3fn+3 


ҚАЛАН (2) 
2 2 
3 3 5 7 
where Хау, Y = 3 1,5 „2, 5 355 ‚4. 
We proceed by seeking an approximate solution of the form 
m+t-1 
yw)= Ý yx (3) 
j=0 
m+t-1 
y= У jay xi t= Лару) (4) 
j=1 


where m and t are the number of collocation and interpolation points used in the method. 
Specifically for this method m = 9,t = 1 and the degree of the polynomial is m + t — 1. 
Equation (3) is interpolated at x = x, and (4) is collocated at 


x= Ха+т: v= 0,2 , 1,5,2,5, 3,2 3 
and 4 which leads to the following non linear system of equations of the form 
mt+t-1 
Р(х) = > ej х} 
j=0 
mt+t-1 т+е-1 
Ро) = У Лер =h У By Of -Боу) (5) 
j=1 j=0 


When using Maple 17 (Mathematical software) to determine the unknown parameters ос; 


and В, in (5), we obtain the following. 


Xo= Yn 


18027 71277 з 2732 4 6989 , 1369 6 161 ,, 


bo =| 775608: 72268042: 108072: 1 5400044: 324005: ' 189012 
ДЕ 18 + Е P 
7560h7 85058 
jc 8200 , 8192(8658) , 2048(1745) , 8192(1316) „ , 4096(575) „ 
1 405405һ 1216215h2 57915h3 289575h^ 173745h5 
16384(73) , 4096(35) „ 65536 
4054056 4054056 3648645h8 
7560 , 28494. 42783 , 33713, 7605 2(989), 69 


md = __ = 6 — 
б = | 1801: 1270027 360m. 825042 72705! + 315h6 ^ ` 50076 





9 


8 
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2052), 4(10338) , 1(16867), 4(1425) 2(6805) , 16(463) , 
9455 


&-| 135h 405h2 135h3 6754 405h5 


г AOD og CA P| 
135һ6 1215A8 
, | 1(1890) , 1(16137), 1(13785),  1(24463) „ 1(6115) „ 1(867) 
B 60h 180h2 120h3 300h* 180h5 105һ6 
1(65) 8 

= 18 P 
60h5 t 135h8 








ІШ 




















ee с” „ _ 8(6606) , |2(1659), 8(1522), 4(789), 32(58), 463), 
2 | 315h 945h? 45һ3 225h4 135h5 315h6 ` ` 315Л6 
128 i 
2835h8 
_ EC 2 1011178) з 120033), 1(18821); 1(5017), | 20761) ,, 
3 324h 48612 648h3 810h4 4865 567h5 








_ (61) e, 16 J 
162h9 ` ` 729h8 








p= 4..4(4842) 4. 1(1257) 1 41202) | 2(655) „ 16(51) „ 12(59) ja 
2 1155h 3465h? 165h3 825h* 495h5 1155h® 11556 
64 
та 
10395h8 


К | 1(4215) „ 1(8541) , 1(3921), 1(15203)  1(4215) „  1(671) „ 
14040h 1404082 4680h3 2340014 1404012 8190Л6 
1(57) 4 
-— 18 + v) 
4680h5 5265h8 








(6) 
where l = (x — Xn) 
у(х) =X Yn + h |Bofa + Baf „з Pf + Baf „үз + Война + Bsf „5 + Bsfn+3 
4 4 2 2 2 2 
+ Baf а + Pf] (7) 
Equation (6) is substituted in equation (7) to obtain our continuous formula. Also 
evaluating (7) at x = хл, v = 2 1,5 ,2, = ‚3 E and 4 to yield eight integrators to form 
our hybrid block methods as follows: 
LE Уң 
10000000 0000000 1 
010 0 0 о o o=] |о о о о о Ы 
001000 о o| |o 0 0 о о о 0 1/3 
0001000 Off¥Yn+z}_}0 0 0 O 0 O о 1|| Yn 
000 01 0 0 02.42 000 00 0 0 112,3 
00000100), | 0 00000 0 ay | 
00000010), |Д |0000 00 0 ту 
00000 о 0 17% 000000 о 14) "2 
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985667 


400400 
46741504 


18243225 
4848 


1925 
46235648 


18243225 
1841200 


729729 
63488 


25025 
502544 


200475 
48234496 


18243225 


+h 


6797493 


63078400 
973 


7425 
27081 


246400 
5552 


51975 
29675 


266112 
27 


275 
272629 


950400 
47104 


51975 


Let A = 


әәссос«оссоск 
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1372587 


` 1192755200 


8413 


` 737100 
4203 


` 358400 
4213 


` 368550 
178775 


` 15095808 
99 


` 9100 
64141 


` 4147200 
24242 


184275 


ooocooocoro 
оооооноо 
оооонороо 
оооноооо 


оонооороо 


72842607 


` 22937600 
340769 


` 113400 
244719 


` 89600 
39386 


` 14175 
797975 


` 290304 
3897 


` 1400 
2809513 


` 1036800 
43072 


` 14175 


ске‚осоосоосососо 











CÍ!ooococcococ 

















14005219 89105481 72842607 2085791 
5734400 45875200 22937600 4587520 
102397 72607 16082 91943 
42525 37800 141750 204120 
61343 360477 1881 8327 
22400 179200 1600 17920 
128144 8098 2272 11462 
42525 4725 2025 368550 
641525 271175 50135 1247875 
217728 193536 36288 2612736 
527 2133 306 71 
175 1400 175 ` 280 
2243563 ` 916153 194089 368039 
777600 691200 129600 1866240 
145408 9896 32768 11584 
42525 4725 14175 С 25515 
ae 0 0 16196113 
91750400 
ЕЕ 0 0 59977 
340200 feta 
000 00 0 0548 fn-3 
358400 Ғ 5 
0 0 0 00 ODER ш 
d 85050 || fn-2 
ТТ 615535 Е 
3483648 
247 fn-1 
0 0 0 0 0 = 
1400 hel 
ao: ð 0 0 2202641 үг! 
12441600 
0 0 0 0 0 21958 
6075 
(8) 
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By multiplying (8) by the inverse of AØ and rearrange it in Butcher Table as 




























































































C A 

3 | 16196113 985667 72842607 14005219 89105481 72842607 2085791 6797493 1372587 

16 | 91750400 400400 22937600 5734400 45875200 22937600 4587520 63078400 1192755200 

1 59977 46741504 340769 102397 72607 16082 91943 973 8413 

4 340200 18243225 7113400 42525 7 37800 141750 7 204120 7425 7 737100 

3 63341 4848 244719 61343 360477 1881 8327 27081 4203 

8 358400 1925 7789600 22400 ` 179200 1600 ` 17920 246400 ` 358400 

1 15011 46235648 39386 128144 8089 2272 11462 5552 4213 

2 85050 18243225 14175 42525 4725 2025 25515 51975 368550 

5 615535 1841200 797975 641525 271175 50135 1247875 29675 178775 

8 | 3483648 729729 ` 290304 217728 ` 193536 36288 7 2612736 266112 ` 15095808 

3 247 63488 3897 527 2133 306 71 27 99 

4 1400 25025 71400 175 71400 175 ~ 280 275 79100 

7 | 2202641 502544 2809513 2243563 916153 194089 368039 272629 64141 

8 | 12441600 200475 ~ 1036800 777600 ` 691200 129600 1866240 950400 4147200 

1 1058 48234496 43072 145408 9896 32768 11584 47104 24242 
6075 18243225 14175 42525 4725 14175 ` 25515 51975 184275 

1 1058 48234496 43072 145408 9896 32768 11584 47104 24242 
6075 18243225 14175 42525 4725 14175 ` 25515 51975 184275 


(9) 


The Table (9) satisfies Runge-Kutta conditions for solution of first order ODEs since 


(i) y ау = Ci 
pa 
(it) 2 


The method (8) is formally given as Runge-Kutta type method as 
529 12058624 10768 36352 2474 
k, + k, К, + k, ak 
э енын 12150 18243225 14175 42525 4725 
"o 7" 18192, 2896 ӨШЕДІ; 12121 


SG tay pea ae А Ко ы ЕТ puerum a 8 "EE 9 
14175 25515 51975 368550 
k, = Ті; у,) 











16196113 - 985667 7284260 14005219 89105481 


+ 
" 367001600 ` 1601600 ^ 91750400 ^ 22937600 * 183500800 ` 
á ‚ 3285531 , _ 2085791, 6797493, 13732587 
11468800 ^ 18350080 ’ 252313600 ` 4771020800 





k, = >, 
2 f Xs 16 
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59977 11685376 340769 102397 72607 

k, + k, k, + k, Е 
an 1360800 18243225 453600 170100 151200 
k,=f| 4 и 8041 91943 793 8413 
+ ks k, + 8 9 
28350 816480 29700 2948400 




















63341 1212 244719 61343 360477 
k; + k, k, + k, st 
k= fix 5% PET. 1433600 | 1925 358400 89600 716800 
6 ch um 1881 8327 27081 4203 
— k, ->k + kg — 9 
6400 71600 985600 ` 1433600 


15011 11558912 19693 32036 4049 
k, + k, k, + k, а 
Б = 1 340200 18243225 28350 42525 9450 
з= |х, +=, у, +А 
2 568 5731 1388 4213 
+ К; „+ 3 
2025 51030 51975 1474200 


615535 , , 460300, _ 797975, 641525, 271175, 
-— 13934502 ` 729729 ^ 1161216 ^ 870912 * 774144 ` 
M 50135 1247875 , , 29675 , 178775 

145152 ^ 10450944 ’ 1064448 ° 60383232 ° 


247 15872 3897 527 2133 
k, * k, К, + k, 2 
5600 25025 5600 700 5600 
153 71 27 99 
6 7 + К, ko 
350 1120 1100 36400 




















9 











5 
ks = f Х, tg 














k,=f x + Óh, y, +h 








2202641 125636 2809513 2243563 916153 
k, + k, k, + k, А 
x 49766400 ! 200475 4147200 3110400 2764800 
? , 194089, | 368089 | | 272629 | 64141 
518400 ^ 7464960 ` 3801600 “ 16588800 ` 


k; =f 7 











529 12058624 10768 36352 2474 
1 + k, k, + k, 5 
p 12150 18243225 14175 42525 4725 
o = f| х, +h, y, th 
à 8192 k 2896 11776 12121 


+ + 
14175 ^ 25515 ' 51975 ? 368550 ` 





(10) 
3.0 Consistency and Stability of the block method 
Thus obtain the normalized form of (8),the first characteristics polynomial of the normalized 
matrix will be expressed as 


p(R) = det[RA® — AY] 
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e 
e 


p) = |А 


oo occ oc он 
осососо-то 
оо оооноо 
оо ооноо о 
оо оноо 0о о 
оонооооо 
оноо оороо 
н ооооооо 
огососососо 
O oOo0 00 SS 
O-O 0000 O DS 
Ooo О O 
OO DÜ OO DO SO 
әсоссосос«сос«ос«сосососо 
әсоосос«сос«ос«ссос«с 
к к ҥш кш кш нн 

II 

© 


A'(A—1)20 
A, = 0,22 = 0,23 = 0,4, = 4, А5 = 0,26 = 0,4, = 0, А8 = 1 
From definition 1.0, the newly hybrid block method (8) is zero stable and also consistent 
since the order of all the integrators are uniform order 9 > 1 
4.0 Numerical Experiments 
The following examples are used to confirm the efficiency of our method 
Example 4.1 


1 
у' = 20x? — 20y + 2x, y(0) = 3 h-005 0<х<10 


Exact Solution: у(х) = x? + Зе 20х 


Example 4.2 

у = -у, у(0) = 1 һ-0.045 0<х<10 
Exact Solution: y(x) = e~* 

Example 4.3 (SIR Model) 

The Susceptible Infected Recovery (SIR) model is an epidemiological model that computes 
the theoretical number of people infected with a contagious illness in a closed population 
over time. The name of this class of models derives from the fact that they involve coupled 
equations relating the number of susceptible people S(t), number of people infected /(t), 
and the number of people who have recovered R(t),. This is a good and simple model for 
many infectious diseases including measles, mumps and rubella. It is given by the following 
three coupled equations. 


ds 
ae y= Bis 


dt |. 
Ши I — yI + BIS i 
7 Vi +B (0 
и Rt yl 
where y, y and Б are positive parameters. Define y to be, 

y=S+I+R (ii) 


when solutions in (i) are substituted in (ii) we have 
y = u(1 — y)t then 
y'=u(-y) y@=p (Ші) 
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Taking и = 0.5 and attaching an initial condition y(0) = 0.5 (Юга particular closed 
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population), we obtain 


у' = 0.5(1 — y), 


y(0) 2 0.5 


Exact solution: y(t) = 1 — 0.5e 9:52 


Example 4.4 
у = A(sinx — y), 
Exact Solution: y(x) — 2 


y(0) = 0 


2 . 
sinx — 





2241 


h = 0.05 





2241 
Absolute errors, with stiffness ratio A = 100 for fixed step size h = 0.01,0 € x < 1.0 


—Ax 





cosx + e 


2241 


Table 1: Approximate solution of Example 4.1 at k=4 


0 < x < 1.0 













































































Mesh Exact solution Present method k= 4 | Absolute 
values error 
0.05 0.125126480400000 | 0.125126480300000 | 1.0 x 10-10 
0.10 0.055111761070000 | 0.055111761050000 | 2.0 x 10-11 
0.15 0.039095689460000 | 0.039095689450000 | 1.0 x 10-11 
0.20 0.046105212960000 | 0.046105212980000 | 2.0 x 10-11 
0.25 0.064745982330000 | 0.064745982360000 | 3.0 x 10—11 
0.30 0.090826250730000 | 0.090826250760000 | 1.0 x 10-10 
0.35 0.122803960700000 | 0.122803960700000 |................ 
0.40 0.160111820900000 | 0.160111820900000 |................ 
0.45 0.202541136600000 | 0.202541136300000 | 3.0 x 10-10 
0.50 0.250015133300000 | 0.250015133200000 | 1.0 x 10-10 
0.55 0.302505567200000 | 0.302505567200000 |.................. 
0.60 0.360002048100000 | 0.36000204800000 1.0 x 10-10 
0.65 0.422500753400000 | 0.422500753400000 |................ 
0.70 0.490000277200000 | 0.490000277200000 |................ 
0.75 0.562500102000000 | 0.562500102000000 |................ 
0.80 0.640000037500000 | 0.640000037400000 | 1.0х10-10 
0.85 0.722500013800000 | 0.722500013800000 |................ 
0.90 0.810000005100000 | 0.810000005000000 | 1.0х10-10 
0.95 0.902500001900000 | 0.902500001800000 | 1.0х10-10 
1.00 1.0000000001000000 | 1.0000000000000000 | 1.0 x 10? 








DERIVATION OF NINTH STAGE RUNGE-KUTTA METHOD 


Table 2: Approximate solution of Example 4.2 at k=4 
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Mesh values Method [9] atk=4 | Method [8] at k=4 Present method at k=4 
0.05 1.7509 х10714 | -------- 
0.10 1.456 x 1078 1.3173 x 10714 1.0 x 10-13 
0.15 1.3607 x 10714 1.0 x 10-15 
0.2 7.2078 x 10? 1.2558 x 1071 1.0 x 10-15 
0.25 1.3668 x 10714 1.0 x 10715 
0.3 1.2499 x 107? 7.5660 x 10 1^ 1.0 x 10-15 
0.35 8.1734 x 10-14 1.0 x 10-13 
0.4 3.4300 x 107? 3.1670 x 10-13 1.0 x 10-13 
0.45 3.1299 x 1071? 1.0 x 10715 
0.5 6.6468 x 10^? 2.9542x 10713 1.0 x 10-15 
0.55 2.8179 x 10713 2.0 x 10715 
0.60 2.0233 x 10? 2.6783 x 10713 1.0 x 10-13 
0.65 2.5574 x 10713 1.0 x 10-13 
0.7 5.8373 x 10? 2.3961 x 10773 1.0 x 10-13 
0.75 2.7780 x 1071’ 1.0 x 10715 
0.8 4.5984 x 107? 4.2459 x 10-13 1.0 x 10-15 
0.85 4.1166 x 10-13 1.0 x 10715 
0.90 2.3751x 10? 3.9011 x 10-13 1.0 x 10-15 
0.95 3.77160 x 10713 1.0 x 10-15 
1.00 5.2617 x 10? 3.5322 x 10-13 1.0 x 10715 























Figure 1: Error graph of Example 4.2 at k= 4 


== Method [9] at k=4 


[ii Method [9] at k=4 


== Present Method at К-4 
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Table 3: Approximate solution of Example 4.3 at k=4 
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Mesh values Method [5] atk=4 | Method [2] ағ | Present Method at 
=4 =4 
0.1 ESUIASEQTX) е ee 
0.2 3.946177 E(-12) | 200 EC15) |------- 
0.3 818 932ЕСИП- Ге ||ыш---де 
0.4 S436119EC11) |800Е(15) -------- 
0.5 1.92974 Е(-10) 1.20 Е(-14) 1.0 Е (-15) 
0.6 1.87904 E(-10) 1.60 E (-14) 1.0 Е (-15) 
0.7 1.776835 E(-10) |180Е(-14) 2.0 E (-15) 
0.8 1.724676 Е(-10) | 2.30 E (-14) 1.0 E (-15) 
0.9 1.847545 Е(-10) | 2.40 E (-14) 2.0 E (-15) 
1.0 3.00577 E(-10) 2.90 E (-14) 2.0 E (-15) 








Figure 2: 





Error graph of Example 4.3 


== Method [5] at k=4 


== Method [2] at k=4 


== Present method at k=4 
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Table 4: Approximate solution of Example 4.4 at k=4 










































































Mesh Method [9] at Method [8] atk=4 | Present method k = 4 
values k=4 

0.05 5.9679 x 104 4.2826125 x 10-8 
0.10 1.8070 x 10-2 1.5485 x 10-4 8.27482906 x 10% 
0.15 8.1938 x 10°5 5.581025 x 10-9 
0.20 1.0210 x 10? 7.5021 x 10-35 5.543383 x 10? 
0.25 1.1954 x 10-6 3.7376 x 10-11 
0.30 1.3225 x 10°3 3.8056 x 10-6 2.53 x 10-13 

0.35 1.0512 x 10° 2.00 x 10:15 

0.40 4.3646 x 103 1.2499 x 10-6 1.00 x 10-15 

0.45 6.6184х 107 | --------- 

0.50 7.8860 x 10-4 1.9415 x 10-8 1.00 x 10-15 

0.55 1.0247 х108 | --------------- 

0.60 4.4568 х 10-4 9.3785 x108 ^ | ----------------- 
0.65 1.4945 х108 | ---------------- 
0.70 5.7728 x 10-4 4.7574x10? ^ — |------------- 
0.75 1.3142 x 10? 1.0000 x 10:15 
0.80 1.9057 x 10? 15626x10?  — | ---------------- 
0.85 8.2737 x 10-10 2.0000 x 10:15 
0.90 3.4431 x 104 2.4271x 1010 ^ | ------------- 

0.95 1.2804х1011 | --------------- 

1.0 1.9459 х 104 1.1724 x 10-11 1.0000 x 10-15 

















Figure 3: Error graph of Example 4.4 


== Method [9] at k=4 


== Method [8] at k=4 


== Present Method at k=4 
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5.0 Discussion of Results 
We observed that the Ninth stage Runge-Kutta method performed excellently well with 


the four problems tested with the method. This shows that our method is good and can be 


used to solve accurately any model of the form y' = f(x,y), y(a) = p. (see Tables 1, 
2,3,4 and figures 1,2 and 3) 
6.0 Conclusion 


We want to conclude that the newly block integrator (8) is of uniform order 9, zero stable, 


consistent and self starting and after reformulating into Runge-Kutta type method, the 


results obtained from the method converges more excellently than the existing method. 


(see figures 1, 2 and 3). Although the cost of implementation is high when compare with 
method [2]. 
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